Variational calculation of the limit cycle and its frequency 
in a two-neuron model with delay 
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We consider a model system of two coupled Hopfield neurons, which is described by delay differ- 
ential equations taking into account the finite signal propagation and processing times. When the 
delay exceeds a critical value, a limit cycle emerges via a supercritical Hopf bifurcation. First, we 
calculate its frequency and trajectory perturbatively by applying the Poincare-Lindstedt method. 
Then, the perturbation series are resummed by means of the Shohat expansion in good agreement 
with numerical values. However, with increasing delay, the accuracy of the results from the Shohat 
expansion worsens. We thus apply variational perturbation theory (VPT) to the perturbation ex- 
pansions to obtain more accurate results, which moreover hold even in the limit of large delays. 

PACS numbers: 82.40.Bj, 84.35. +i, 02.30.Ks, 02.30.Mv 
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I. INTRODUCTION 



Feedback in biological systems has received increased 
attention in recent years In particular, the role of de- 
layed recurrent loops in models of population dynamics, 
epidemiology, physiology, immunology, neural networks, 
and cell kinetics has been studied extensively 0] ■ Neural 
network systems are complex and large-scale nonlinear 
dynamical systems, and the dynamics of a delayed net- 
work are yet richer and more complicated Hopfield 
Q proposed a simplified model of a neural network in 
which each neuron is represented by a linear circuit con- 
sisting of a resistor and a capacitor, coupled to other neu- 
rons via nonlinear sigmoidal activation functions. From 
this model, he derived a system of first-order ordinary 
differential equations to describe the network dynamics. 
Extending Hopficld's model, Marcus and Westervelt Q 
considered the effect of including a temporal delay in the 
model to account for finite propagation and signal pro- 
cessing times. 

In networks of real neurons, delays occur at the synap- 
tic level due to transmitter release dynamics and the 
integration time of post-synaptic potentials at the den- 
dritic tree level where post-synaptic potentials have a 
finite conduction speed to the soma, and in the axons 
due to the finite axonal conduction speed of action po- 
tentials 0]. It is well-known that time delaycan cause 
an otherwise stable system to oscillate f?,^J^ and may 
lead to bifurcation scenarios resulting in chaotic dynam- 
ics Uii • On the other hand, delayed feedback permits 
the control of chaos [IH, where it can be used to stabi- 
lize unstable periodic orbits in chaotic attractors [13, 14]. 
Experimentally, time-delayed chaos control was success- 
fully applied, for instance, to electronic oscillators p^ . 
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mechanical pendula [l^|, lasers [l7|, and chemical sys- 
tems [l8| . Furthermore, a recently proposed scheme for 
the treatment of neurological disorders employs delayed 
feedback in order to efficiently desynchronize the activ- 
ity of oscillatory neurons [l^. Therefore, finite delays are 
an essential property of any realistic model of a neuron 
population [20|. 

In the vast majority of cases, information about a phys- 
ical system can only be obtained by means of numerical 
or analytical approximation methods. Numerical meth- 
ods constitute a powerful and effective tool to describe 
even extremely complicated physical scenarios. Never- 
theless, their accuracy is not always superior to that of 
analytical approximations, and usually more insight into 
the physical principles that govern the system is obtained 
by pursuing an analytical approach. Often, perturbation 
expansions are easily accessible, but they are usually di- 
vergent and need resummation. A recently developed, 
powerful method to perform such a resummation is vari- 
ational perturbation theory (VPT), which has been suc- 
cessfully applied in various quantum or statistical field 
theories [H, [H, H, HI . A first apphcation of VPT 
in the field of deterministic nonlinear dynamics is found 
in Ref. [26], while the present work extends the use of 
VPT for the first time to a system described by delay 
differential equations (DDE's). 

In Sec. [ni we introduce the two-neuron model and the 
system of DDE's that we consider. The results of a lin- 
ear stability analysis of the model system are reported in 
Sec. IIIIl and it is shown that a limit cycle emerges via a 
supercritical Hopf bifurcation when the delay exceeds a 
critical value. In Sec. lIVi the Poincare-Lindstedt Method 
is applied to derive the perturbation expansions for the 
delay-induced limit cycle and its angular frequency. In 
Sec. |Vl we apply the Shohat expansion to the perturba- 
tion series of the limit cycle and its angular frequency 
as a first crude resummation approach. In Sec. IVIl we 
resum the perturbation expansions using VPT, which al- 
lows us to improve the quality of our results significantly 
and to obtain results which are reasonable even in the 
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FIG. 1: Numerical solutions of the system of DDE's ((TJ, with ai = —1, a2 — 2 and r*^' = r'^' = r. For this choice of 
parameters, the critical value of the delay is ro = 7r/4 ~ 0.7854 ... . In a) the delay is r = 0.7, and the origin is a stable fixed 
point. In b) the delay exceeds the critical value: r — 0.8. In this case, the origin is unstable and the trajectory approaches a 
limit cycle. In both cases the initial conditions are ui{t) — 0.2, U2{t) = for t G [— r, 0]. 



limit of large delays. 



Setting 



II. MODEL 

Neural circuits composed of two or three neurons form 
the basic feedback mechanisms involved in the regulation 
of neural activity [l^l . Many researchers have used bifur- 
cation analysis and numerical simulations in order to ana- 
lyze a system of two Hopfield-like neurons with discrete or 
distributed time delays [IJ H H [s^, IM 112, HI, H HI ■ 
In this investigation, we apply analytical approximation 
methods to a two-neuron system with delay, described 
by the coupled first-order DDE's 



dui{t) 
dt 

dU2{t) 

dt 



-|-aitanh[u2(i-T(2))] (1) 
~U2(i)+a2tanh[ui(t-r(i))]. (2) 



Here, ui and U2 denote the voltages of the Hopfield neu- 
rons and T^^^ and t*-^-* are the signal propagation or pro- 
cessing time delays, while ai and a2 describe the cou- 
plings between the two neurons. 



III. LINEAR STABILITY ANALYSIS 

The system of DDE's ([1]), ^ has a trivial stationary 
point at ui = M2 = and we first analyze its stability. 
Near the equilibrium point, linearizing the DDE system 
yields 



ui{t) = -~ui{t) + aiU2{t - T^"^^) , 

U2{t) = -U2{t) + a2Ui{t - T^^^) . 



(3) 
(4) 



u(i) 



(5) 



in the last equation, where A is a complex number, and 
ci and C2 are constants, we get a nontrivial solution if 
and only if 



aia2e 



-A(r< 



0. 



This equation has been analyzed in detail in Ref. [2^ 
For aia2 < —2 the conditions of Theorem 2 in Ref. [26 
are met. Defining t = (t*^^^ + r(^^)/2 and 



1 

2loq 



sm 



2tj| 



aia2 



2j^ 



J = 0, 1, 2, 



(7) 



where cjq = \/|aia2| ^ 1, this theorem states that: 

• If T e [0, To), then the zero solution of ([T]), ([2]) is 

asymptotically stable. 

• If r > To, then the zero solution of (P), dH) is un- 
stable. 



• Tj, with J = 0, 1,2 
of©, @. 



. , are Hopf bifurcation values 



Furthermore, Theorem 3 in Ref. [28| states that the Hopf 
bifurcation at t = To is supercritical. Note that iujQ is 
the solution to ([6|) when t = tq, and the period of the 
limit cycle at the Hopf bifurcation is thus Tq = In jujQ. 



IV. POINCARE-LINDSTEDT METHOD 

Figure [T] shows numerical solutions of the system of 
DDE's ID), ^ for the two cases in which the delay t is 
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either smaller or greater than its critical value. Below 
the critical value tq of the delay r no periodic solution 
exists, while above t — tq there is such a solution. We 
now consider the case rf^^^ = t*-^-* = r, aia2 < —2 and 
seek to calculate the period and trajectory of the pe- 
riodic solution approximatively. To this end, we apply 
the Poincare-Lindstedt method [36[. Since a supercriti- 
cal Hopf bifurcation occurs at r = tq, we assume that 
the amplitude and frequency of the new periodic states 
are analytic in e = \/t — tq and expand them as 



-^ + ^vl^\^-u.oro). (19) 

UJQ UJQ 



Imposing the initial conditions V^^\o) = Aq, V'2^'^''(0) — 
Bq on the periodic solution V*^°^ (^) , we find the general 
solution to the system of homogeneous differential equa- 
tions mi), (HH) as 

Vl°\0 = ^0 cose + Soaisin(c^oT-o) sin ^, (20) 



rio). 



u(<) = eV{t) ^ e [u(°)(t) + eU(i)(<) 



(8) 
(9) 



V^'\0 = i?ocose 



Ao 



ai sin(a;oTo 



■ sin e . 



(21) 



It is convenient to rescale the argument of these functions 
so that they become periodic with period 27r. We thus 
introduce the new independent variable ^ according to 



and we write 



U(t) = V(0 . 



(10) 



(11) 



The periodic solution V(e) to ([12]), uM can only be de- 
termined up to an arbitrary phase. Without loss of gen- 
erality we can thus choose = in (HH, which 
fixes the phase of the zeroth order solution, at least up 
to a shift of TT. 

In general, to order e", we have to solve the system of 
differential equations 



Applying the perturbation expansion ([5]) to the system 
of DDE's (H]), ([2]) and performing the change of variables 
(fTU)) . PT|) . we obtain 



dvt\0 



(22) 



UJQ 



UJQ 



dvm 

di 
dV2{0 



d^ 
in which 

a(e) 



ai 



— tanh{ey2K-a(e)]} , (12) 



^tanh{eyi[e-a(e)]} ,(13) 



(23) 



uj{e)T = w(e)(ro + e^) 

uJoTo + eujiTo + e^(wo + ^2To) + 



(14) 



The delayed variable V^i/2[e — Q;(e)] is written as 
V[e-a(e)]=V(")(e,a)+6VW(e, 



where the inhoniogeneity f '^"^ (e) is determined by the 
solutions to previous orders. Since we require that the 
solution V*^") (f) be periodic in ^ with period 2tt, we can 
impose certain conditions on the inhomogeneity f {^) . 
Namely, we demand that f ^^"^ not contain terms that 
would lead to non-periodic solutions for V^")(e), i.e., 
f(")(^) must not contain secular terms. In order to iden- 
tify the conditions that must be satisfied by f^"^(e), we 
expand V'"' (C) and f {£_) as Fourier series: 



) + ..., (15) /l//"^(e)\_ 



corresponding to the expansion in ([S]) , which is equivalent 
to 

V(0=V(")(e)+eV(i)(e) + ... . (16) 

In order to take into account (fT4|) . each term in the ex- 
pansion of V(e — a) is expanded as a Taylor series: 



V2^"\0 



/l*"^(0 

ft\0 



E 

k=l 



E 

fc=i 




a 



(") 
i,fe 

(n) 
^2,k 



COS fc^ 



cos fc^ -I- 




sin 



(24) 



Pl,k 
gin) 



sin 



(25) 



V(^')(C-cooTo) 



(17) 



dv(j')(e') 



+ 



By inserting the expansions ([M|) . (|25p into the systems 
of equations (P^ . we find that the coefficient of the 
the terms with fc = 1 in the inhomogeneity f ^^"^ (e) must 
satisfy the conditions 



Applying the expansions for V(e) and V(e — a) to (fT2|) . 
()13|) . we obtain to zeroth order in e 



a2sm{LjQTo)a\{ + fr^.i = 



dVl°\0 
d^ 



^ + ^Fi")(e-c.oro), (18) 



aa'i - 0.2 sin(woTo)/3i'i = 



(26) 
(27) 



UJQ 



UJQ 



The derivation of these two conditions is demonstrated 
in the appendix. 
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n 




2 


4 

2 + 7r 


4 


4(341 + lOSvr) 
27(2 + 7r)3 


6 


8(73843 + 4077371 + 58327r^) 
729(2 + 7r)5 


8 


1440729464 + 37r(359606308 + 928145677r + 83980807r^) 
98415(2 + 7r)7 


10 


2(1885638326848 + 97r(193375795408 + 37r(22966214893 + 47r(952738307 + 629856007r)))) 

13286025(2 + n)'^ 


12 


(48294520193761504 + 37r(17432699637100336 + 37r(2577825095210584 
+7r(596088219927028 + 7295935409444l7r + 38093690880007r^))))/(8370195750(2 + tt)") 


14 


-(137083613818976067424 + 37r(56352224911533618320 + 37r (9835626348748269040 + 37r( 949 130678879606440 
+37r(54285350368574420 + 7r(5287281140608997 + 2285621452800007r))))))/(1129976426250(2 +7r)^^) 


16 


(290578164278923471719089408 + 97r(44452665928743252091582336 + 37r(8868376426577693217600640 
+37r(1013305929108995195501920 + 97r(24272564564656648301080 + 7r(3331148075811324207916 
+7r(270489187825118497343 + 99835945058304000007r)))))))/(111054083171850000(2 +7r)^^) 



TABLE L Expansion coefficients for the angular frequency of the hmit cycle for ai = —1, a2 = 2 up to order e 



After this general result, we now consider the first- 
order expansion of the system (fT2|) . (fT3|) . Taking into 
account the result (PO)) . (PT|) and the choice Bq = 0, we 
obtain the inhoniogeneity f '^^^ to be given by 



/(')(0-Aoa>i(ToCose ^ 



To 



sin^ 



(28) 



and 



-Aouji 



02(1 + To) 



UJQ 



+ 



To 



oi sin(a'oTo 



sin(a;oTo) cos^ 



■sin^ 



(29) 



Thus, according to the conditions ([26]) . ([27]) . we must 
demand 



2Aoa;iTo 
ai sin^ (woTo) 



and 



AAoLUijl + To) 
ai sin(2(jJoTo) 



0.(30) 



We must thus have either Aq = or = 0. If we 
choose ^0 = Oi we only obtain the trivial solution. Thus, 
we choose uji — 0, and the coefficient Aq is yet to be 
determined. The solution for V*^^^ (^) is then simply given 
by the solution to the homogeneous system: 



A\ cos^, 
^1 



■sinC, 



ai sin ( Wo To) 

where A\ is to be determined in higher orders. 



(31) 
(32) 



Expanding (fT^ . up to order while taking into 
account the zeroth- and first-order result, we obtain the 
inhomogeneity f (^) as given by the expansion (pS)) . For 
the first component we have 



(2) 
^1,1 



(2) 



Pl,3 



4afwo 



f ^0(^0 + T0W2) 

^o(l+To)t^2 , 



12a? 



Wo 



^0 



And for the second component we have 

ylo(l + To)w2 

4 



J2) 



(2) 
"2,3 



-02 sin(woTo) 
-02 sin ( Wo To) 



Wo 

—2- -I- ^o(wo + T0W2) 
4wo 



-4o 



^3 

-a2sin(woTo)^ [2cos(2woTo) - 1] , 

^3 

^02 sin(woTo)-—- 2- [2cos(2woTo) + 1] 

IzWq 



(33) 
(34) 
(35) 
(36) 

,(37) 
(38) 

(39) 
(40) 



while all other coefficients vanish. Imposing the condi- 
tions (|26p . (P?]) on the inhomogeneity f^^-'(^), we obtain 
the system of equations 

A^(l -l-ai -l-w^) - 8a?wo(wo -f W2T0) = 0,(41) 
Agwo(H-a?-|-w^) + 8a?wo-l-8a?(l-|-ro)w2 = 0.(42) 
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(n) 


A: = 1 


= 3 




k = 1 


A: = 3 


n = 


4 

A/3(2 + 7r) 





n = 








n = 2 






n = 2 







81(2 + 7r)-V2 


27(2 + 7r)3/2 


27(2 + 7r)3/2 


(n) 


it = 1 


A; = 3 




k = 1 


A = 3 


n = 








n = 


4\/2 
A/3(2 + 7r) 





n = 2 





10^6 


n = 2 


^(436 + 937r) 




27(2 + 7r)3/2 


81(2 + 7r)5/2 


27(2 + 7r)3/2 



TABLE 11: Fourier expansion coefficients of the limit cycle for ai = —1, 02 = 2 up to the second order in e. 



Its solutions read 



UJ2 



Ao 



■ To 



2 ' 



(1 + al + Lu^){l + To + luItq) 



(43) 
(44) 



Choosing the sign of Aq to be positive fixes the phase 
of our zeroth order solution definitively. This procedure 
can easily be carried to higher orders, where only even 
orders lead to nonzero terms for both the corrections to 
the angular frequency ujn and the expansion of the limit 
cycle V(")(f). Expanding ^,^to order e^", we find 
the coefficient ^2(71-1) a-nd the correction LU2n from the 
conditions (HH), ((27)) . 

From here on, we consider the choice of parameters 
fli = —1, 02 = 2. These parameter values lead to luq = 1, 
To — tt/4: and the solution (|^^ reduces to 



UJ2 = 



v/3(2 + n) 



(45) 



Table U shows the first eight nonvanishing corrections to 
the angular frequency. Note that the signs of the ex- 
pansion coefficients ojn alternate and that their absolute 
value grows rapidly. This indicates that the perturba- 
tion series for w is a divergent Borel series. Table HIl 
shows the expansion coefficients of the first two nonvan- 
ishing orders of the Fourier expansion of the limit cycle 
as given by ([M)) . Figure [5] a) shows the first eight orders 
of the perturbatively calculated angular frequency uj^'^\ 



N 
ri=0 



(46) 



as a function of e. Note that odd and even perturba- 
tion orders yield results which are respectively smaller 
and larger than the numerical values. For small values 



of the delay, the perturbative solution is in good agree- 
ment with the numerical data. However, as e grows, the 
perturbative solution becomes unacceptable. Figure[2]b) 
shows an example of the perturbatively calculated limit 
cycle given by 



it) 



N-1 



(47) 



where we count the order N of our perturbation expan- 
sion such that in the A^th order we obtain the A'^th nonva- 
nishing corrections lu2N and V(2(^-i))(^). For the value 
e = 1 chosen in Fig.[2]b), the limit cycle can still be ob- 
tained with good precision from the perturbation series 
(|T7)) and as in the case of the angular frequency we ob- 
serve that the perturbative approximations approach the 
numerical result in an alternating manner. However, as e 
increases, the perturbative solution ([T7)) becomes useless 
as in the case of the angular frequency. Thus, if we want 
to obtain analytical results for larger values of the tem- 
poral delay t, we must resum our perturbation series. In 
the next section, we apply a Shohat transformation to 
the perturbative results for both the angular frequency 
cj(^) and the limit cycle u'-^'>{t). 



V. SHOHAT EXPANSION 

Now, we resum our perturbative results by performing 
a Shohat expansion. This method was first introduced 
for calculating the period of a Van der Pol oscillator in 
Ref. [131 and it has been conjectured that the expansion 
yields results which are valid for all values of the per- 
turbation parameter [37i, i38i] . Furthermore, it has been 
stated that the Shohat expansion is succesful when the 
periodic solution to the differential equation in question is 
of softening type, i.e., the angular frequency lo decreases 
with e |39i] . which is the case for our system as is evident 
from Fig. [5] a) . 
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ui{t) 



FIG. 2: Perturbative results for the angular frequency u) and the limit cycle {ui{t), U2{t)}. In a) the angular frequency is 
shown as a function of e. The dashed lines represent the perturbative results as given by (|46|) and Tab. |T] Numerical results are 
shown by dots. In b) the limit cycle {ui{t), U2{t)} is shown for e = 1. Dashed lines represent perturbative results according to 
(|47|) : the numerical result is shown by the solid line. 



The basic idea of the Shohat expansion is to map the 
perturbation parameter e € [0, oo) to a new parameter 
/i G [0, 1). In order to perform the resummation of the 
angular frequency, we introduce the new expansion pa- 
rameter /i according to the transformation suggested by 
Shohat in Ref. 37] and thus set 



1 



(48) 



where we explicitly take into account that the pertur- 
bation series for the angular frequency depends only on 
even powers of e. Inverting we have 



(49) 



1 — /i 



We now obtain the Shohat expansion of our perturbative 
result by replacing in ([^ according to the last identity 
and re-expanding the series in /i up to order /i^. The 
Shohat expansion of the angular frequency is thus given 
by 



N 



Ul. 



{N} 



ri=0 



fe=0 



1 



(50) 



The resummation of the limit cycle (j47p is performed in 
a similar manner and we obtain 



AN) 



N-l n 
n=0 k=0 



ri — 1 

k 



^V(2("-'=))(i). (51) 



Finally, in order to evaluate the resummed angular fre- 
quency and limit cycle for a certain value of e, we replace 
fi in ([50)1 and ([?T|) according to 



Figure [3] a) shows the angular frequency after Shohat 
resummation as a function of the delay parameter e. We 
find that, if we go to sufficiently high orders, the re- 
summed expansion yields reasonably good results for all 
values of e. Figure [3] b) shows an example of the limit 
cycle after resummation. Note that for the value of the 
delay parameter in Fig.[3]b) the perturbative result prior 
to resummation would be completely useless. 

Figure 0] a) shows the convergence of the angular fre- 
quency obtained from the Shohat expansion versus the 
perturbation order N for different values of the tempo- 
ral delay. For small values of the delay, the convergence 
seems to be exponentially fast, at least up to the eighth 
order. For larger delays, the convergence appears to be 
less regular. In order to examine the convergence of the 
limit cycle results, we consider the error measure 



<5( 



N) 



/;°+^dt ||u(0-uW(0| 



To+T 



dt II u( 



(52) 



where we rescale the argument of our analytic solution so 
that its period is identical to the period of the numerical 
solution and shift the phase of the analytic solution ac- 
cording to the phase of the numerical solution. Figure [4] 
b) shows the convergence of the results for the limit cy- 
cle. As in the case of the angular frequency, the results 
from the Shohat expansion and their convergence with 
the perturbation order are best as long as the delay is 
not too large. In the next section, we thus use a more ef- 
ficient method to resum the perturbation series. It yields 
accurate results already in low orders, allows us to ob- 
tain more precise results, and its convergence depends 
less crucially on the size of the delay parameter. 
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TV = 2 7V = 4i; 
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ui{t) 



1.5 



FIG. 3: Angular frequency and limit cycle after Shohat resummation. In a) the angular frequency is shown as a function of 
e. Results from the Shohat expansion as given by (|50|) are shown by dashed lines. Numerical results are shown by dots. The 
inset shows a magnification of the interval 4 < e < 5. In b) the limit cycle is shown for e = 2. Dashed lines represent results 
from the Shohat expansion as given by (|5ip; the numerical result is shown by the solid line. 



VI. VARIATIONAL PERTURBATION THEORY 

In this section, we improve the resummation of the per- 
turbation series of the angular frequency and the limit cy- 
cle by applying VPT to the perturbation series (|46)) and 
(j47|) . This method is based on a variational approach due 
to Feynman and Kleinert [2l|, which has been system- 
atically extended to the nonperturbative approximation 
scheme now called VPT [2l,MlMl2i. 



A. Basic Principles 

VPT is capable of converting divergent weak-coupling 
into convergent strong-coupling expansions and has been 
applied successfully in various fields, such as quan- 
tum mechanics, quantum statistics, condensed matter 
physics, and the theory of critical phenomena. In fact, 
the most accurate critical exponents come from this the- 
ory [ioj, as verified by recent satellite experiments [4l| . 
First applications of VPT in the field of Mark ov p rocesses 
and nonlinear dynamics are found in Refs. [47l . |4^ and 
Ref. [ll|j respectively. 

The convergence of VPT has been analyzed up to very 
high orders for the ground-state energy of the anhar- 
monic oscillator 



1 



V{x) = ^ u^x^ + gx" 



(53) 



and was found to be exponentially fast [42, Ej]. This 
surprising result has been confirmed later by studying 
other physical systems and was proven to hold in general 
[HiHJ]. Furthermore, the exponential convergence seems 
to be uniform with respect to other system parameters. 
The variational resummation of perturbation series thus 



yields approximations which are generically reasonable 
for all temperatures ^44 45 ] , space and time coordinates 
0, l47l l48l, m agnetic field strengths [49] , coupling con- 
stants [2^, [50, |5l| J spatial dimensions [5^], etc. 

VPT permits the evaluation of a divergent series of the 
form 



N 



(54) 



and yields a strong-coupling expansion of the generic 
form 



M 



fid) - 9"^' E ^" 



9 



-2m/q 



(55) 



m=0 



Here, p and q are real growth parameters and character- 
ize the strong-coupling behavior. Introducing a scaling 
parameter k, which is afterwards set to one, Eq. (|54[) can 
be rewritten as 



N 



n=0 



(56) 



Applying Kleinert's square-root trick [23[, i.e. setting 

K = K^/T^^ (57) 

with 



(58) 



in (j56p . the variational parameter K is introduced into 
the perturbation series: 



N 



f^''H9) = Ea„g"ifP-"'? (1 +gr)(''-"' 



9)72 



K=l 



(59) 



8 




N N 

FIG. 4: Convergence of the angular frequency and the hmit cycle after Shohat resummation. In a) the logarithm of the relative 
deviation of the angular frequency as given by (|50|l from the numerical value and in b) the logarithm of the error measure for 
the limit cycle as given by (|52|l are shown versus the perturbation order. In both a) and b) different symbols indicate different 
values of e (dots: e — 1.6; squares: e — 2.0; triangles: e = 3.0; diamonds: e = 4.0, upside-down triangles: e — 5.0). 



The Taylor series of (1 + gr)" vi^ith a= {p — nq)/2 reads 



N-r 



{1+grr 



; = 1 



k=0 



1 



0(g 



N-n+l\ 



(60) 



where the generalized binomial coefficient is defined by 

r(a + i) 



T{k + l)T{a - k + 1) 



(61) 



The series (l60l) is truncated after k 



N 



n since the 

N 



original function Z*-^-* (g) is only known up to order g 
As a result of this truncation, the function /^^-'(g) be- 
comes dependent on the variational parameter K: 



n=0 k=0 ^ ^ ^ 



(62) 

k 



The influence of the variational parameter is then opti- 
mized according to the principle of minimal sensitivity 
[sst . i.e., one evaluates the function ((62|) at that value 
of the variational parameter K for which it has an ex- 
tremum or turning point. In the following, we set g — 
in dini) and (|T7)l . 



B. Resummation of the Angular Frequency 



when the physical quantity in question has a strong- 
coupling expansion of the form ((55)) [23l [2^ . There- 
fore, we first consider our numerical data for the angular 
frequency in the case of large delays and determine the 
growth parameter p and q in (j55p . To this end, we analyze 
our numerical data in two steps. First, in Fig. [5] a), we 
plot our numerical results for Incj versus Ing = ln(T — tq). 
Fitting our data to a function of the form 



/(ln.g) = p/qlng + Infeo , 



(63) 



we find p/q — —0.9997 and 6o = 1.565. We expect the 
growth parameters to be integers and thus set p/q = — 1. 
For large delays, the leading asymptotic behavior of u) is 
thus given by w ~ g"^. In order to determine not only 
the ratio of p to g but the individual values of the growth 
parameters, we then fit our data for gu) to a function of 
the form 



/(5"') = &o + M 



-2/q 



(64) 



which is shown in Fig. Ob). The numerical results from 
the fit are: bo = 1.571, bi = -2.7, and q = 1.993. Thus, 
we assume q = 2 and from our previous result we then 
have p = —2. In order to determine 6o and bi numerically 
with better accuracy, we can now perform a hierarchy of 
approximations to order M by fitting guo to functions of 
the form 



fig- 



M 



(65) 



We now apply VPT to obtain an improved resumma- 
tion of the angular frequency VPT is applicable 



From this procedure we obtain the more precise numer- 
ical values 6o = 1.57081 and bi = —2.66. Now, we can 
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a) 



In w 




In^ g-'^/m-^ 

FIG. 5: Angular frequency for large delays [y/r — to £ [50, 100]). In a) the logarithm of the angular frequency is shown versus 
the logarithm of the delay parameter g. Numerical data are represented by dots; the solid line represents a fit of the data to 
a function of the form (|63p . In b) the product of the delay parameter and the angular frequency is shown versus the inverse 
square of the delay parameter. Numerical data are represented by dots; the solid line represents a fit of the data to a function 
of the form ((64)) . 



introduce the variational parameter K to the perturba- 
tion series (|46l) according to ([S^ with p = —2, q = 2: 



(66) 



n=0 fc=0 ^ ^ ^ 

To first order we obtain 

(1) f (2 + 7r)(2j^^-l)-4g 



K^{2 + tt) 



which has a minimum at 



x(i) = Wl 



4.9 

2 + 71 



(68) 



Evaluating (|67p at the optimized value of the variational 
parameter then yields 



2 + 7r 



45 + 2 + TT 



(69) 



In the limit of large delays, g ~* 00, we thus have 

4^)^(5, if «)^4V^ + 6iV^ (70) 



with 



and 



foW = ^ « 1.28540 



= -1.6522. 



16 



(71) 



(72) 



To second order, Eq. (|66p yields 



1 



27ii:6 

108 (2 
+ 9 TT 



27 {3K^ - SK"^ + 1) 



3K' 



9 



,4(341 



(73) 
1087r)" 



(2 + 7r)3 



Since this has no real extremum in the variational pa- 
rameter K, we look for roots of the second derivative. 

In general, in order to optimize the influence of the 
variational parameter, we first look for minima or max- 
ima oi uj^r^{g , K), and if those do not exist, for positive 
roots of higher derivatives. In each order iV, the op- 
timized variational parameter i^T^^-* is thus determined 
from the condition 



or 



dK 

<pJ^%{g,K) 



<PK 



= 0, ., 



(74) 



In cases where a certain derivative has several positive 
roots, we choose the one which is closest to the optimized 
value from the previous order if(^-i). The iVth order 
VPT approximation of the angular frequency is then ob- 
tained by evaluating (j66p for the value of the optimized 
variational parameter: 



^vpt(5) 



(75) 



Returning to ([75)l . we find that for 



< 3(2 + .)[24+12. + 5735(2T^] ^ ,, ,,, 
^ - 587 - 1447r ^ ' 
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FIG. 6: Angular frequency and limit cycle from VPT. In a) the angular frequency as given by (|75p is shown as a function of 
e for the orders A'^ = 1, 2, 8 of VPT (orders three through seven would lie very close to the curve for A'^ = 8). Dots represent 
numerical values. The inset shows a magnification of the interval 4.8 < e < 5. In b) the limit cycle is shown for e = 2. Dashed 
lines represent the results from VPT as given by (|9ip . The numerical result is shown by the solid line. 



ujypr^{g, K) has two positive turning points: 

~ (2) _ ^60 + 157r(4 + tt) + 60(2 + 7r)g ± 2ri 



with the abbreviation 



3(2 + 7r) 



(77) 



r/ = V2T^ (78) 
X v/9(2 + 7r)3 + 72(2 + nYg - (587 - 1447r)52 . 

Comparing ([77)1 to we find that K_' is closer to 

if^i) and thus evaluate (El) for K = K^^^ to obtain 



(2) , £>(2)s 



27(2 + 7r)2 



(79) 



[15(2 + 7r)(2 + 7r + 4.g) -2?]]^ 
X I4772 - 42(2 + 7r)(2 + TT + 4g)r] + (2 + tt) 

X [117(2 + nf + 936(2 + irfg + 4(1061 + mn)g^]^ . 

However, for delay parameters exceeding the value of g 

~ ("2") 

given in (|76p . we cannot use Kj_ since in this case 77 be- 
comes imaginary. Thus, if we want to consider the limit 

of large delays, we must optimize the variational param- 

f 2) 

eter by considering the third derivative of ujyprj,{g, K), 
which turns out to have two positive roots for all posi- 
tive g: 



(2) _ ^180 + 457r(4 + tt) + 180(2 + n)g ± p 



with the abbreviation 



3^/2(2 + tt) 



, (80) 



P = V2 -h TT 



(81) 



v/513(2 7r)3 -t- 4104(2 + irfg + 16(5137r - 724).^ 



Again, K_ ' is closer to the first-order solution, and we 
set ^ ^(2)^ ^^^^.^^ 



54(2 + 7r)2 



[45(2 + 7r)(2 + 7r + 4g) -p]3 
|p2 - 72(2 + 7r)(2 TT + 45)p + (2 + 7r)[1323(2 + -nf 



+ 10584(2 + 7r)2g + 16(2771 + 13237r)52]| . 
Expanding the last result in g^^, we obtain 



with 



Ug — 



27(2-H7r)3 



90 + 457r - v'(2 + 7r)(5137r- 724) 



2047 + 18367r - 72a/(2 + 7r)(5137r - 724) 



and 



l(2) 



243(2 -Ftt) 



90 + 457r - v/(2 + 7r)(5137r- 724) 
437 



(82) 



(83) 



(84) 



1.23174 



(85) 



x<'63213W — '^-^^ — + 1627r 
[ V 5137r - 724 



2 + 7r 



5137r- 724 



82 



-I- 426v/(2 + 7r)(5137r- 724) - 16193 \ w -1.1229 . 



It thus turns out that the second order approximation 
for the leading and subleading large-delay coefficient is 
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a) 



-2 

-3 
-4 
-5 
-6 
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-8 



b) 



In- 



In (5 



(W) 
VPT 



4 5 6 7 

N 



5 

iV 



FIG. 7: Convergence of the angular frequency and the hmit cycle after resummation with VPT. In a) the logarithm of the 
relative deviation of the angular frequency as given by (|75p from the numerical values and in b) the logarithm of the error 
measure for the limit cycle as given by (I52|l are shown versus the perturbation order. In both a) and b) different symbols 
indicate different values of e (dots: e = 1.6; squares: e — 2.0; triangles: e — 3.0; diamonds: e = 4.0, upside-down triangles: 
e = 5.0). 



actually worse than the first order one. However, the 
results in higher orders turn out to be improved approx- 
imations. For fixed values of the coupling constant, the 
procedure in higher orders is analogous to the first and 
second order, where the roots of the first, second, or third 
derivative of ujyplj,{g,K) have to be determined numer- 
ically. Furthermore, in order to obtain the coefficients 



6 Wand bi^'K 



we expand the derivatives of u!yprj,{g, K) 



\TL g and the variational parameter K as 



(86) 



in order to carry out the optimization procedure. 

Fig. [6] a) shows our VPT results for the angular fre- 
quency versus the delay parameter e. The first order 
result is already in good agreement with the numerical 
results for a wide range of delays and is far superior to 
the first order result from the Shohat expansion (com- 
pare Fig. [3] a)). Figure [7] a) shows the convergence of 
our VPT results for five different values of the delay. 
The accuracy of our VPT results improves with increas- 
ing order; however, not as regularly as in the case of the 
Shohat expansion for small delays. Figure [5] a) shows 
a comparison of the eighth order results obtained from 
Shohat resummation and VPT. In particular, for larger 
values of the delay, the results from VPT are far superior 
to the ones from Shohat resummation. Table IIIII shows 
our results for the leading large-delay coefficients 6o and 
the subleading coefhcient &i; again, the convergence is 
not montonic, but we do observe a general trend towards 
improved results in higher orders. 



C. Resummation of the Limit Cycle 

We now proceed to perform a variational resumma- 
tion of the limit cycle following the approach of Ref. [5^ . 
To this end, we consider the perturbation series of each 
coefficient in the Fourier expansion of V(^) as given by 



(87) 



N-l 

AN) _ (2«) „ 

1/2, fc ~ "i/2,fcy ' 

n=0 
N~l 

1/2, fc ~ "i/2,fcy • 

n=0 



We introduce the variational parameter K into the per- 
turbation series for A^y2f. and i?|^^j,in the same way as 
for the angular frequency, and obtain by applying 
to the Fourier expansions ([87]) 



A 



{N) 

1/2, fc, VPT 
N-l 



(89) 



n=0 



and 



B 



(N) 



1/2, fc, VPT 
N-l 



k=0 



n=0 fe=0 ^ ^ ^ 



(90) 



Instead of optimizing ([59)1 and ([5D|) according to the prin- 
ciple of minimal sensitivity, we obtain our VPT result for 
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9 10 



FIG. 8: Comparison of the eighth order results for a) the angular frequency and b) the limit cycle obtained from the Shohat 
expansion and VPT. The relative deviations of the analytical results from the corresponding numerical values are shown versus 
the delay parameter (Shohat expansion: squares; VPT: circles). 



the limit cycle more easily by evaluating all Fourier ex- 
pansion coefficients for that value of the variational pa- 
rameter K which was determined through the optimiza- 
tion procedure of the frequency, i.e., our VPT result for 
the limit cycle reads: 



^2,VPT V^-' 



E 

k=l L 



A. 



(N) 

2,fc,VPT 



sin 



cos fc^ 



, (91) 



''2,fc,VPT 

where iC'-^"^-' is determined from the condition (|74p and 
we use ^(^-1) instead of K^^\ since the A^th term in 
the series for V(^) is a correction of order g^^^ . 

As an example, we consider the lowest order in which 
we can perform the VPT resummation of the limit cycle. 
To order g our solution for V(5) reads 



4cos^ 
v/3(2 + ^) 
2^/3 



■5\/3(116 + 337r)cosC 



81(2 + 7r)5/2 
[cos3C-7sin3^]l +0(5^), 



27(2 -t- 7r)3/2 ' 
4^/2sinC J \/6(436 -f 937r) sin^ 
V3(2T^"-''l 81(2 + 7r)5/2 (93) 
2^/6 



[5 cos 3^- sin 3C] \ + 0{g 



27(2 + 7r)3/2 



Introducing the variational parameter K according to 



89]), (|90l), we obtain 



1 ) cos ^ g 
if4 



/s:V3(2 + 7r) 
5\/3(116-f 337r) ^ , 2^3 



81(2 + 7r)5/2 



■ cos ^ ^ 



27(2-F7r)3/2 



(94) 

cos 3^—7 sin 3^] 



v.. 



(2) 



2, VPT 



4(2^^2-l)^/2sinC g 



2^/6 



r\/6(436-f937r) . ^ 
\ 81(2 + 7r)5/2 27(2 + ^)3/2 



[5 cos 3^ — sin 3^] 



The optimal value of the variational parameter for the 
angular frequency to first order is given by (j68p . Insert- 
ing this value into ([M)) . ([55)1 . we find the following VPT 
result for the limit cycle: 



(2) 



^(0 = 



1 



;{l08(2-f 7r)2cosC 



27^3(2 + 71) (2 + 7r + 4g)2 
g [(1148 + 6997r) cos^ - 6(2 + 7r)(cos 3^ - 7 sin 3^] } 



(2) 



<0 = - 



(96) 

;{l08(2-f 7r)2sin$ 



2.VPT vs; - 27^6(2 + 7r)(2 + n + Agf 
+ 5[(1292 + 77l7r)sin^-f 6(2-K7r)(5cos3C-sin3^)] } . 

(97) 

The procedure in higher order is analogous. Figure [5] b) 
shows our VPT results for the limit cycle for e = 2 up 
to the eighth order. Figure [7] b) shows the logarithm of 
the error measure (|52p for the VPT limit cycle versus 
the order TV for different values of e. In Fig. [5] b) the 
accuracy of the eighth order results from the Shohat 
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iV 


1 


2 


3 


4 


5 


6 


7 


8 


numerical 




1.2854 


1.23174 


1.56495 


1.59507 


1.61990 


1.61806 


1.61139 


1.54478 


1.57081 


6'^) 


-1.65 


-1.12 


-2.72 


-2.79 


-3.05 


-3.03 


-2.98 


-2.21 


-2.66 



TABLE III: Leading and subleading coefficients for the large-delay behavior of the angular frequency. 



expansion and VPT are compared. Again, we find that 
our VPT results are more reliable than those from the 
Shohat expansion, especially for larger delays. 



VII. SUMMARY 

We have performed a perturbative calculation of the 
limit cycle and its frequency in a two-neuron model with 
delay. A Shohat resummation of the respective pertur- 
bation expansions yields results which are in good agree- 
ment with numerical values but whose accuracy decreases 
drastically with larger values of the delay parameter. Re- 
summing the perturbation series with VPT yields more 
uniformly converging results, which are reliable even in 
low orders, and furthermore permits the extraction of 



the leading large-delay behavior with sufficient accuracy. 
The present work constitutes the first application of VPT 
to a system of DDE's. Moreover, it establishes a method 
for the variational resummation of perturbatively calcu- 
lated limit cycles in nonlinear dynamical systems. 
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APPENDIX A: ELIMINATION OF SECULAR TERMS 



We now demonstrate how the conditions (P7|) are obtained by considering the Fourier decompositions of the 

periodic solution and the inhomgeneity. Inserting ([24]) and ((25|) into the system of equations ((22|) , ([23]) and comparing 
coefficients of sin and cos fci^ in both components, we obtain the following system of four equations: 

(Al) 
(A2) 
(A3) 
(A4) 



(n) 
«l.fc 


+ - 


(n) 


cos(fcwoTo) -t- 


An) 


sin(fcwoTo) 


(n) 


= 0, 


UJQ 


Wo 


Wo 


"l,k 

Wo 


- ka'rl - 


An) 

Wo 


cos(/cwoTo) — 


(") 
aia2,fc 

Wo 


sin(fcwoTo) 


Pl,k 


= 0, 


(n) 

Wo 


+ kb^l ' 


(«) 

Wo 


COs(fcwoTo) + 


Wo 


sin(fcwoTo) 


(n) 
'^2,fe 


= 0, 


"2,k 
Wo 


'^"'2,k 


dn) 
a2b\l 

Wo 


cos(fcwoTo) — 


(n) 

Wo 


sin(fcwoTo) 




= 0. 



It turns out that for A: > 1 the coefficients a^"\ b^"' can be uniquely determined for any inhomogeneity, i.e., for 
arbitrary Q![,"\ /Si"''- For k > 1 the solution of the system (|Aip - (jA4p is 

'^i"fe ^ ~ ^^oM"2)(^o + k^^l) ~ aiuol sin(fcwoTo)(2fca2"^ - (1 + fc^)wo/32"fc ) + aiwo cos(fcwoTo) (A5) 

X (2ag - 2fcwo/32 + (1 - k^W^a^^l) + (wo + wg) [sin(2fcwoTo)(/3g - fcwoag) + cos(2A:woro)(ag + fccJo/^S)] } , 

^i"fc = ^^oal"fc)(^^o + k^^l) - fli^o sin(fcwoTo)(2fc/3^"2 + (1 + fc^)^oa2"fc) + °^i^o cos(fcwoTo) (A6) 

X (2/3("2 + 2A;woa^";^ + (1 - k^Hn'^^l) + (c^o + w^) [cos(2A;woro)(/?g - fcwoag) - sin(2A:woro)(ag + fc^o/^S)] } , 
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-ai 



(A7) 



^(") _ 



cos(fcwoTo)(2ag - 2fcwo/?S + " fc')'-'oaS) + sin(2fcwoTo)(/32 - fc^o^S) + cos(2fcc^oTo)(a2 + fc^o/^^J)] } ' 
= TtI (/32'2 + A:c^o41)(^o + + (c^o + ^o') [— sin(fcc^oTo)(2fc/3|",^ + c^o(l + fc')al"fc) (A8) 

cos(fcwoTo)(2/3i;'2 + 2fcwoal"fe + (1 - fc')woMr2) " sin(2fccjoro)(4?fc + fc^o/3^"2) + cos(2fcworo)(/3^"2 " fc^oa2?fc)] } > 



where 



D = 2 + 2tj^(l + k^) + L0^{1 + /c^) + (wo + t^o) [2(1 - fc^^o) cos(2woTo)/wo - 4fc sin(2fcToWo)] 



(A9) 



Note that D vanishes for k = 1. We must thus reconsider the system (|Aip - (|A4| for the case k — 1 and it turns 

out that OLi^\ f3[^^ must satisfy certain conditions for a solution to exist. For fc = 1, we add 02 sin(Li;o'n3) x ()Alll to 
^4p and subtract ujq x (jA3p from a2Cos(woTo) x (jA2p . Using the identities aia2 = — (t^o + 1) and wq — cot{u}oTo), 



we obtain the two conditions (|26p . (jST)) that must be satisfied by the inhomgeneity f'^"'(^). Imposing (|26|) . ((27| on 
f(")(^), we obtain the foUowing solution to the system of equations (|Aip - (IA4p for fc = 1: 



b'f'} = cos(tJoTo) 



a2 cos{ujQTo)a\ { + 



(n) 

Q2a 
a2 



a2sin(u;oT-o) ' 



621 = a2sm(woro) { - { sm(Ci;oTo) cos(woTo) + "2,1 cos (woTo) . 



(AlO) 
(All) 



Here, the coefficients a^"i, aj"^ are undetermined and follow from the initial conditions. We set aj^"/ = An and 



(«) 



,(") 

*2,1 



0. 
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